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ABSTRACT 

The distance-redshift relation determined by means of gravitational waves in 
the clumpy universe is simulated numerically by taking into account the effects 
of gravitational lensing. It is assumed that all of the matter in the universe 

c^ ! takes the form of randomly distributed point masses, each of which has the 

identical mass M^. Calculations are carried out in two extreme cases: A ^ 
k> \ GMl/(? and A ^ GMi/c^, where A denotes the wavelength of gravitational 

H \ waves. In the first case, the distance-redshift relation for the fully homogeneous 

and isotropic universe is reproduced with a small distance dispersion, whereas 
in the second case, the distance dispersion is larger. This result suggests that 
we might obtain information about the typical mass of lens objects through the 
distance-redshift relation gleaned through observation of gravitational waves of 
various wavelengths. In this paper, we show how to set limitations on the mass 
Ml through the observation of gravitational waves in the clumpy universe model 
described above. 

Subject headings: gravitational lensing - gravitational waves - dark matter - 
distance scale 
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Introduction 



Owing to recent developments of observational techniques, it is expected that detection 
of gravitational waves will be achieved in the near future. Laser interferometer gravitational 
wave detectors, such as the Tokyo Advanced Medium-Scale Antenna (TAMA-300), the Laser 
Interferometer Gravitational- Wave Observatory (LIGO), Variability of Irradiance and Grav- 
ity Oscillations (VIRGO) and GEO-600, are in operation, and other projects, such as the 
Large-scale Cryogenic Gravitational wave Telescope (LCGT), LIG02, the Laser Interfer- 
ometer Space A ntenna (LISA), t he Decihertz Interferometer Gravitational Wave Observa- 



tory {DECIGO; ISeto et al.l (1200 ll )). and the Big Bang Observer (BBO), are in the planning 
stage. In preparation for future breakthroughs associated with these projects, it is useful 
to discuss the information provided by gravitational wave data. Here we focus on inhomo- 
geneities of our universe that can be investigated with gravitational wave observations. 

Many cosmological observations suggest that our universe is globally homogeneous and 
isotropic, in other words, that the time evolution of the global aspect is well approximated 
by the Friedmann-Lemaitre(FL) cosmological model. However, our universe is locally in- 
homogeneous. The effects of these inhomogeneities on the evolution of our universe or on 
cosmological observations draw our attention as one of the most important issues in cosmol- 
ogy. However, it is difficult to directly observe the inhomogeneities with optical observations 
because most of the matter in our universe is not luminous. One useful means of obtaining 
information about aspects of inhomogeneities is to examine gravitational lensing. Gravita- 
tional and electromagnetic waves are subject to gravitational lensing effects caused by inho- 
mogeneously distributed matter around their trajectories. Therefore, we may find features 
of inhomogeneities through both observational and theoretical studies of the gravitational 
lensing effects in astronomically significant situations. 

There are various candidates for the dark matter such as weakly interacting massive 
particles (WIMPs) and massive compact halo objects(MACHOs). Several observations have 
been undertaken to constrain the mass density of compact objects, f^co; direct searches for 
MACHOs in the Milky Way have been performed by t he MACHO and ER OS collaborations 



through microlensing surveys. The MACHO group (jAlcock et al.l |2000| ) group concluded 



that the most likely halo fraction in form of compact objects with masses in the range 0.1- 



^E-mail:c_ni_yoo@sci. osaka-cu.ac.jp 
^E-mail:knakao@sci. osaka-cu.ac.jp 
■^E-mail:furusaki@sci. osaka-cu.ac.jp 
'*E-mail:Takahasi@th. nao.ac.jp 



-3- 



IMq is of about 20%. The EROS team concluded that the objects in t he mass range from 
2 X 10~'^Mq to IMq cannot contribute more than 25% of the total halo (jAfonso et al.ll2003l ). 
However, the universal fraction of macroscopic dark matter could be significantly different 
from these local estimates. Millilensing tests of ga mma-ray bursts deriv e a limit on flco of 
^co < 0.1 in the mass range from 10^ to 10^ Mq ( iNemiroff et al.l 120011 ). Multiple imaging 
searches in compact radio sources derive the limit on fico as fico i$ 0.013 in the mass range 
from 10^ to 10^ Mq ( IWilkinson et al. Il200ll ). The mass range of these observational tests 
are limited, and thus other methods to investigate the mass fraction of macroscopic compact 
objects in wider mass range are needed. 



In this paper, we show that it is possible to extract information about the properties of 
macroscopic compact objects from the observational data of gravitational waves by analyzing 
the gravitational lensing effects due to these compact objects. For this purpose, we consider 
an idealized model of the inhomogeneous universe which has the following properties. First, 
this model is a globally FL universe. Second, all of the matter takes the form of point 
masses, each of which has the identical mass Ml. Finally, the point masses are uniformly 
distributed. 

One of the main targets of gravitational- wave astronomy is binary compact objects. 
Their gravitational waves have much longer wavelengths (A) than the optical, and further- 
more, these are coherent. In typical optical observations, the wavelength (A ~ 1/im) might 
be much shorter than the Schwarzschild radii of lens objectsO Thus, we usually analyze 
the gravitational lensing effects on electromagnetic waves by using geometrical optics. In 
contrast, we need wave optics for the gravitational lensing of gravitational waves, since the 
wavelength of gravitational waves may be comparable to or longer than the Schwarzschild 
radii of lens objects. For example, the wave-band of LISA is typically 10^^ — 10^^ cm; a point 
mass of 10^ — IO^Mq has a Schwarzschild radius almost equal to the wavelengths within 
that wave-band. When we consider gravitational lensin g effects on the g ravitational waves 
with GMl/c^ < A, we have to take the wave effect (JNakamural Il998l ) into account. In 
fact, remarkable differe nces between the extrerne case s A ^ GMl/c^ and A ^ GMl/c^ have 
already been reported (JTakahashi fc Nakamurall2003l ). 



In order to obtain information about uniformly distributed compact objects by investi- 
gating the gravitational lensing effects in observational data, we focus on the relation between 
the distance from the observer to the source of the gravitational waves and its redshift. The 



■''It does not conflict with any observational result that there might be lens objects of Schwarzschild radii 
much smaller than 1/zm. However, for example, supernovae which might be the smallest and yet very bright 
optical sources, are larger than the Einstein radii of such small-mass lens objects, and thus the lensing effects 
due to those will be negligible. 
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distance is determined using the information c ontained in the ampl itudes of the gravita- 
tional waves from the so-called standard sirens (JHolz &: Hughes! l2005l ). Alth ough there are 
sever al similar analyses using the optical observation of Type la supernovae (JHolz fc Linder 
20051 ). gravitational lensing effects on gravitational waves can add new information about 
the properties of the lensing compact objects to that obtained by optical observations, by 
virtue of the wave effect, through which we can gain an understanding of the typical mass 
and number density. 

Here we assume that the redshif t of each source is in dependently given by the observation 
of the electromagnetic counterpart (JKocsis et al.ll2006l ) or the waveform of the gravitational 
waves (JMarkovidll993l ) . The luminosity distance d\ fr om the obser ver to the source is given by 
the waveform of the gravitational waves, as follows (ISchutzlll986l ). Since the orbit of binary 



will be quick ly circularized by a gravitational radiation reaction (jPeters fc Mathews! Il963 



Peters! 1 1964! ). the effect of ellipticity is negligible when the emitted gravitational radiation 
becomes so strong that we can observe it. Then, the two wave modes emitted from a binary 
of two compact objects with masses mi and m2 are given by 



h. 



h. 



2Mtr^{ 



di 
7r/)2/3sin(27r/t) 



d] 



(cos^0 + l). 



2cos^, 



(1) 



(2) 



where / and 6 are the frequency and the angle between the angular momentum vector of the 
binary and the line-of-sight, respectively, and A^ chirp is the redshifted chirp mass, defined by 

{17111712)^^^ 



Alchirp = (1 + z)Mchirp = (I + Z) 



(mi + m2)^/^ 



Assuming a circular orbit, we have 

dt 



967r^/3 



/"/^-MJ^L- 



(3) 



(4) 



Eliminating A^chirp and 9 from equations ([I]), (El), and (jlj), we can obtain di if there is no 
gravitational lensing effect on the gravitational waves. However, if the gravitational waves 
are gravitationally lensed, their wave forms are changed. The gravitational lensing effects 
may lead to an incorrect estimation of di, and thus we have to know how the gravitational 
lensing changes the waveforms of the gravitational waves. 

This paper is organized as follows. In ^ the basic theory of gravitational lensing is 
introduced. We discuss the lensing probability in ^ Then we show the calculation method 
and results for the distance-redshift relation in §§1] and [51 Finally, ^ is devoted to the 
conclusions and summary. Throughout the paper, we adopt the unit of G = c = 1. 
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2. Review of gravitational lensing 

In this section, we introduce basic equations for gravitational lensing. 

2.1. Geometrical optics 

In the case of a poiri t-mass lens with mass Ml, the bending angle vector cx{$,) is given 
by JSchneider et al.lll992h 

1^1 



a(0 = 4Mi^, (5) 



where $, is the impact vector (see Fig. [T]). For convenience, we consider the "straight" line A 
from a source to the observer and define the intersection point of A with the lens plane as 
the origin of the lens position C, and the ray position 7. We note that, in general, the vector 
^ is not orthogonal to the line A. However, if the geometrically thi n lens approximation i s 



valid, we can regard the vectors ^, C,, and 7 as parallel to each other (JSchneider et al.lll992l ). 



We define Ds, Dl, and Dls as the angular diameter distances from the observer to the 
source, from the observer to the lens, and from the lens to the source, respectively. Using 
simple trigonometry (see Fig. [1]), we find the relation between the source position 77 and ^, 



^=^^-^L5«(0- (6) 



D 



L 



Equation (jH]) becomes 



= ^7-^L5«(7-C) (7) 



D 



L 



in terms of ^ and 7. 

Here we introduce the dimensionless vectors defined by 









X = 


4 


and I = 


77 

^0' 


where ^0 


is 


the Einstein radius 


given by 

^0 = 








h'^^s 


LS 


and 















(8) 



(9) 



(10) 
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Obscrvcr 



Fig. 1. — Geometry of gravitational lensing by a point-mass lens. The impact vector ^ 
represents the relative position of a light ray on the lens plane to the lens position, and ex. is 
a vector whose norm is equal to the bending angle. 



From Figure [T], we find the relation 






(11) 



The source will have two images at position {x^,x^). Substituting equation ([5]) into equa- 
tion ([6]) we find 

x± = ]^[Vf^Tl±l), (12) 

where x± = \x±\ and / = \l\. The total magnification factor /i is given by 



/i = l/i+l + |^_| 



det 



dl 
dx: 



det 



dl 
dx_ 



f + 2 



(13) 



2.2. Wave optics in gravitational lensing due to point mass lens 

For simplicity, let us consider gravitational waves that propagate in asymptotically fiat 
spacetime with the metric 



ds' 



[1 + 2U)dt^ + (1 - 2f/)dr2, 



(14) 



where U is the gravitational potential and we assume |f/| ^ 1. If the wavelength of the 
gravitational waves is much smaller than the typical curvature radius of the background 
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spacetime, the equation for the gravitational waves is equivalent to that of scal ar waves. 
Assuming monochromatic waves of an angular frequency cj, we have (jPeterslll974j ) 



(V^ +ij^)(j) = 4ij^U(j). 



(15) 



For convenience, we introduce the amplification factor F defined by 

0H' 



(16) 



where is the plane wave with no lensing effects, and (p^ is the wave that is undergoing 
lensing effects. Here note that 0^ is no longer a plane wave, although it was initially. In the 
case of the point-mass lens with mass M^, the gravitational potential is given by U = —M^/r. 
Then we have the solution of equation (fT5|) in the form 



F{uj; y) = exp 



nw w 
X + ^2^" 



r 1 



-w 



iFi 



-w,l;-wy' 



(17) 



where F and iFi are th e gamnia fun ction and the confluent hypergeometric function, respec- 
tively, and w := AuM^ (jPeterslll974l ). Since |Vf/| is much larger than the Hubble parameter 
of the universe, by replacing u with ^^(l + zl), we can use equation (TT71) in cosmological 
situations. 



2.3. The Dyer-Roeder distance 

In the clumpy universe model, the local geometry is inhomogeneous, but its global aspect 
might be well described by the FL universe whose metric is given by 



ds' 



-dT^ + a^{T) 



dR^ 



1 + KR^ 



+ R^dn- 



(1^ 



where K = 1, 0, and —1, and dQ'^ is the round metric. However, the global aspect of the 
clumpy universe is a vague notion since we do not have a definite mathematical prescription 
to identify the global aspect of the locally inhomogeneous universe with the FL universe 
model. Therefore, the above FL universe is a fictitious background universe; hereafter, we 
will refer to this fictitious geometry (universe) as the "global geometry" (or universe). We 
assume that the global geometry is well approximated by that of the FL universe filled with 
dust and with a cosmological constant. 

In the clumpy universe model, the distances Dl, Ds, and Dt,s in equa t ion (El) or equa- 
tion ([7]) are replaced by the Dyer-Roeder (DR) distances (JKantowskil Il969l : iDyer &: Roeder 



19731 ). which are the observed angular diameter distances if the gravitational waves prop- 
agate far from the po int masses. The DR distance with smoothness parameter a = 
( jSchneider et al.lll992l ) in the global universe defined above is given by 



Ddr{Zi,Z2) 



l+Zi 



dz 



:i + ^)' vano(i + zy - fiKo(i + zy + Qao'' 



(19) 



where zi and Z2 are the redshift of the observer and source, and fimo, Hq, and ^ao are the 
present values of the total density parameter, Hubble parameter and normalized cosmological 
constant in the global universe, respectively, and ^ko = ^mo + ^ao ~ 1- Hereafter, the 
subscript means the pr esent va l ue. S o me further discussions of t he propertie s of th e DR 
distance are presen ted in iLinderl (119881 ). ISeitz fc Schneiden (119941 ) , iKantowskil (119981 ). and 
Sereno et all J200lh . 



3. The distribution of point masses 

In order for the clumpy universe to be the same as a homogeneous and isotropic universe 
from a global perspective, the point masses must be distributed uniformly. However, since 
there are no isometries of homogeneity and isotropy in the clumpy universe, the notion of 
"uniform distribution" cannot be introduced in a strict sense. Thus, we assume a uniform 
distribution of the point masses with respect to the global geometry represented by equa- 
tion (I18p in a manner consistent with the mass density of this global universe. This universe 
model is called " the on-aver a ge Fri edmann universe." The assumptions in this universe model 
are discussed in ISeitz et al.l (11994 ) . 



The comoving number density p^ of the point masses M^ is given by 

3"m0-f^0 



Pn 



m7 



StcMt 



(20) 



where p is the average mass density and we have set Oq = 1. We consider a past light cone in 
the global universe. The vertex of this light cone is at the observer, at i? = and z = 0, and 
is parametrized by the redshift z of its null geodesic generator. Then the comoving volume 
of a spherical shell bounded by R{z) and R = R{z + Az) on this light cone is given by 



AV = AttR^^Az. 

dz 

Therefore, the number of point masses in this shell is given by 



AA^ = pnAy 



OiirnO^n ^odii 



2Mr 



R'^Az, 

dz 



(21) 



(22) 
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where 



dR 

dz 




l + Hl^K,{l 



zYD%{z) 



(23) 



where Df{z) is the angular diameter distance from the source of the redshift z to the observer 
in the global universe, which is expressed as 

1 smXiz) 



where 



Dpiz) 



Hoi^ + z) vn^ 

1 sinhXf^;) 



I Ho{i + z) v/=n 



for fixo > 
for r^Ko = 
for Qxo < 



KO 



(24) 



Xiz) 



V\^ 



KOl 



dz' 



Yiz) 



^n,^o{i + z'f- 

dz' 



nKo(i + z')^ + n 



AO 



(25) 



(26) 



/o Vn„,o(l + Z')^ - fiK0(l + Z'y + ^Ao ' 

The average number of point masses in the region y < C/^o < y + ^y of this shell is given 
by 



Piy,z)Ay 



'^<lyAy 



AN. 



47ra2i?2 
Substituting equations (|22|) and (|23|) into equation fl27|) and using equation ([9]), we have 



(27) 



p(y, z)Ay = 3ifo^mo(l + ^) 



1 



H^QKoil 



zfDUz) 



fi„,o(l + Zf - fiA'o(l + ^)2 + ^AO 

Ddr{^,z)Ddr{z,zs) 



-yAyAz, 



(28) 



DoRiO, zs) 

where zs is the redshift of the source and we have replaced each of distances in equation ([9]) 
with the DR distance. We find that p{y, z) does not depend on Ml, and therefore the lensing 
probability for a given source does not depend on Ml. 

We note that the assumption of the distribution of point masses given in equation fl28l) 
does not necessarily guarantee equality between the magnification in the global universe 
defined by 

'Dor{0,z) " 



fXp 



DAz) 



(29) 



and the av erage magnificat i on in the clumpy universe. Seyral previous works have discusse d 
this issue (|Weinberglll976l : IeUIs et allliJ : lciaudell[2nnnl : iRoselbonil : kibble fc Li^boosl ). 
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4. The long-wavelength case, A ^ Ml 



In the case of A ^ M^, the amplificatio n factor F{u}{l + ZL)',y) is almost equal to 



unity (see Fig. [2]) because of diffraction effects (JNakamuralllQQSl ). and we have 



^^^^i^ = \F{u;{l + zl); y) - 1| « 1. (30) 

Thus we assume that a wave (j) can be treated as a spherical wave whose center is at the 
source during all stages of the propagation. Let us consider the situation in which N point 
masses are located at (z^, |/i), 2 = 1,..., N „ where i is the number assigned to each point 
mass in order of distance from the observer. If our assumption is valid, the wave is changed 
whenever it propagates near each point mass as 



N 



'N(, 



where 



i=l 



F, = F{uj{l+Zi)-y,). (32) 

Thus the observed amplitude |0obs| is given by 

N 

I0obs| = n 1^*1101- (33) 

Then, the total amplification |Ftotai| is given by 

N 
^ total 



ni^^i- (34) 



1=1 



We randomly distribute point masses so that the distribution of those is consistent with 
equation fl28|) . First, we divide the spherical region z < zs into A^ concentric spherical 
shells, each of which is bounded by two spheres, z = zi — /\z/2 and z = Zi + /\z/2. We 
take into account only the nearest lens in each shell. As described below, this procedure is 
valid for large numbers A^, i.e., small Az. We find from equation (l28ll that in the ith shell 
Zi — Az/2 < z < Zi + /S.z/2, there is one point mass within the region y <Yi on average, 
where Yi is defined by 

Yi 

p{y,Zi)dy = l. (35) 



Then we randomly put a point mass within the region y < Yi in the ith shell. Here it is 
worth noting that Yi is an increasing function of A^, so by setting A^ large, we can take into 
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5 

y 



X/Ml=10 
VMl=100 



10 



Fig. 2. — Amplification factor F{uj{1 + ZL);y) given in eq. ( ITTj) with replacement uj -^ 
(1 + -Zi^)^^. We plot |F — 1| as a function of y for three cases of \/Ml = 10, 10^, and 10^. 
The redshift of the source zs and of the lens zl are fixed for z^ = 2 and zl = 1. 



account lenses far from the line A. Hence, in the limit of A^ ^ oo, we can all lenses take 
into account. The results do not depend on the value of N, if it is large enough (see Figs. [3] 
and HI). In our calculations, we set A^ sufficiently large. 



Calculating the total amplification factor -Ftotab the observed distance is given by 

DDBiO,Zs) 



a 



obs 



(36) 



total 



We have calculated 2000 samples in two clumpy universe models, the globally Einstein-de 
Sitter (EdS) and globally FL universe with ^ao = 0.7 and i^mo = 0.3. The results are 
depicted in Figures [5] and [6l where we have randomly chosen the sources of the redshift s in 
the range < 25 < 2. 

The distance dispersion tends to vanish as \/Ml increases, and -Dobs is almost equal to 
Df{zs) for sufficiently large X/Ml (see Figs. [7] and [8]). As described in the next section, 
in the case of \/Ml <^ 1, the dispersion is much larger than for \/Ml ^ 1. It can be seen 
from Figures [7] and [S] that the distance dispersion does not sensitively depend on the value 
of fimO) but become large as z increase. 
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Fig. 3. — Observed distance Dohs, depicted as a function of A^ in the case of X/Ml = 10^. 
The sohd hne represents the DR distance in the EdS universe. The dashed hne represents 
the angluar diameter distance in the global universe. 
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Fig. 4. — Same as in Fig. [3l but here the global geometry is given by the FL universe with 
Qxo = 0.7 and Q^q = 0.3. 
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Fig. 5. — Distance- redshift relation of 2000 samples, depicted in the long- wavelength cases 
A/Mi > 1. We have set A^ = 1, 000+ [10, 000 x zs\ in the cases of A/M^ = 10^ and 10^ and 
A^ = 1, 000+ [80, 000 X zs\ in the case of A/M^ = 10^. We have chosen the EdS universe, i.e., 
the FL model with ^ao = and fimo = 1, as a global universe. The solid line represents the 
DR distance in the EdS universe. The dashed line which is hidden behind the data points 
represents the angular diameter distance in the global universe. 
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Fig. 6. — Same as in Fig. [5l but here the global geometry is given by the FL universe with 
^Ao = 0.7 and fimo = 0.3. 
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Fig. 7. — Relation between 6Dohs/ (-Dobs) and X/Ml, where SDohs and (-Dobs) are the rms 
deviation and the average distance given by 100 samples of Dobs in each X/Ml and zs- 
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Fig. 8. — Same as in Fig. [TJ but here the global geometry is given by the FL universe with 
^Ao = 0.7 and ^7^0 = 0.3. 
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5. The short-wavelength case, A ^ Ml 

In the case of A ^ Ml, the value of the amplification factor F(ti;(l + ZL)',y) is much 
larger than unity for small values of y. Thus, we cannot use the linear approximation used 
in the previous section. Instead, here we adopt the geometrical optics approximation. We 
comment that it is difficult to resolve multiple images in the sky using gravitational wave 
detectors, because their angular resolution is ~ l'-l°, much larger than the typical image 
separation. In geometrical optics, lensing does not change the wave form except for its 
amplitude (due to the magnification) and arrival time (due to the time delay) in each image. 
The observed waveform that is given by the superposition of waves from each image has 
undergone amplification and phase shift. In this paper, we assume stationary configurations, 
and the phase shift is not an observational variable. Therefore, we can observe only the 
amplitude of the superposed waves, and we cannot identify the signal as a superposed one. 



5.1. The method of calculation 



In order to calculate the amplification factor, we use the multiple lens-plane method (jSchneider et al 



19921 ). We consider the "straight" line A from a source to the observer, and we put A^ lens 
planes between them so that the straight line A intersects these lens planes vertically. The 
results do not depend on the value of N if it is large enough; we set N = 10^ in our calcula- 
tions. 

In the same manner as in §U we determine the redshifts Zi and randomly situate a 
point mass Ml within the region y < Yi on each lens plane Sj. In the case of A S> Ml, we 
have considered the lensing effects due to all point masses that we situate within the region 
y < Yi. By contrast, in the case of A <^ Ml, we can neglect the lensing effects due to the 
point masses that are far from the ray. Thus we take into account lensing effects due to 
the point masses in the region y < i/max, where y^ax is some positive constant smaller than 
Yi. In other word, we take into account only a few lens planes in which a point mass is put 
within the region y < ^max- Hereafter we denote the number of lens planes that we take 
into account as A^^^:- The validity of this procedure will be supported by the discussion of 
methods for determining an appropriate value of j/max in §5.21 

We denote the lens planes from the observer to the source as Si, S2, . . . , Sat^. For 
convenience, we name the source plane Eat^+i, which is also orthogonal to the line A. We 
define the intersection point of A with Ej as the origin of the lens position ^^ and the ray 
position 7j. Suppose lens positions Cj (^ = 1, • • • , Nl) are given. Then once the ray position 
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7i is given on the first lens plane Si, the ray position 7- is recursively given by 









(37) 



on the jth lens plane S^ (j < A'"^, + 1), where Di = D£)r{0,zi), Dj = D^ji^O, Zj), and 
Dij = Ddr^zi, Zj). For later convenience, we specify the ray positions and the lens positions 
by the following dimensionless quantities: 



Ui 



Then equation (!37|) becomes 



Ui 



Di 



ui- ^ (3ijOL{ui-q^), 



where denoting Diat^+i and Djat^+i by Ds and Dis, respectively, 

D^.Ds 



and 



A. := 



OL{Ui - qj : = 



D,D,s 
Dis 



D, 



«(7^-C^)• 



(38) 

(39) 
(40) 



(41) 



(42) 



Since more than one image appears in the present case, the amplification factor Ftotai 
is determined by both the interference of the rays corresponding to those images and the 
magnification of each ray due to the gravitational focusing effects. We specify a ray from 
the source to the observer by the ray position ui = V(^p^ on the first lens plane Si, where we 
have introduced the label p to distinguish between the rays. The magnification yU(p) for the 
ray p is defined by 

^(^^ = d^o;;^' ^''^ 

where A(p) is the Jacobian matrix of the lens mapping. 



Au 



du 



N+l 



iP) 



dui 



Ut=V 



i-t^(p) 



N 



I 



E^' 



'i,N+l 



dcx{ui - Qi) dui 



i=l 



dui 



dui 



(44) 



Ui=V 



(p) 
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Fig. 9. — Definitions of tlie vectors ckj, C,^ and 7^. 



Tlie amplification factor -Ftotai is then given by 

^total = ^ |/i(p)|^/^exp(i6'(p)) 



(45) 



wliere tlie phase 6(j,) is determined by the time delay T(p) relative to th e ray that does not 
receive any lensing effect, and by the number of caustics n(p) on ray p (JLandau fc Lifshitz 
19621 ). in the manner 

(46) 



TT 



'■= ^T(p) - -n(p). 



The time delay T(p) is given by (ISchneider et al.lll992l ) 



i=i ^ *'*+^ ' 



2-Di,i+i 
and the number of caustics n^p) is give by 



(47) 



1ir=f(p) 



j=l 



(48) 
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where nj(p) is the number of caustics between Sj and Sj+i on ray p and its value is given by 



n. 



i(p) 



for det 



1 for det 



2 for det 



dui 



i+l 



dui 
dui 
dui 



> and tr 



<0 



> and tr 



dui 



i+l 



dui 



>0 



(49) 



dui 



i+l 



dui 



The square of the absolute value of -Ftotai is given by 



total 






<0 



^(g)) • 



(50) 



The first term on the right-hand side of equation ( 15011 represents the she ar and focusing effect 
on th e ray bundle, and the second term represents the interference effect ( JNakamura &: Deguchi 
1999h . 



Finally, we explain the method for finding rays from the source to the observer. If the 
lens positions C,^ {i = 1, ..., Nl) are given, the dimensionless ray position ■u^v^+i on the source 
plane S^r^+i is completely determined by the ray position Ui on the first lens plane Si; 



**JVl+1 = UNl+1[Ui 

Therefore, 'U(p) is a root of the equation 



0. 



(51) 



(52) 



We solve equation (|52|) using the Newton- Raphson method, which requires initial estimates 
for the desired roots (JRauchlll99ll ). In order to find appropriate initial estimates, we consider 
a square region on the first lens plane Si whose center agrees with the intersection point 
of the straight line A with Si (see Fig. ^ and with sides of length 2yinaxCoi) where ^oi = 
\/4:MlDiDis/Ds- In this square region, we place evenly spaced n x n grid points whose 
intervals are equal to ^y^ax^oi/n. Substituting the position vector of each grid point to ui 
on the right hand side of equation flSTl) . we obtain the corresponding ray position xtAf^+i on 
the source plane. Then, we adopt the grid points which lead to 



77Vi,+l 



DsUn 



L + l 



as the initial estimates of the Newton- Raphson method. 



< 1 



(53) 
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In practice, the interval between the grid points must be finite for the numerical calcu- 
lations, so we may fail to find initial estimates that lead to some solutions i7(p). However, the 
magnification factors |/^(p)| of such solutions are so small that the contributions to |Ftotai| 
will be negligible. In order to see this, let us consider the inverse map of equation ( ISTIl . The 
inverse image of S on the first lens plane consists of more than one connected region since 
multiple images of the source appear due to the lensing effects. If a grid point is contained 
in one of the connected regions, the grid point is found as an initial estimate by imposing 
the above criterion given by (1531) . If there is no grid point in a connected region due to its 
smallness, we cannot find initial estimates for the solution V(^p^ included in the connected 
region and cannot find the solution 'U(p) itself. Here we should note that the magnification 
factor associated with a root obtained from initial estimates contained in a connected region 
is almost proportional to the area of the connected region divided by the area of S. There- 
fore the magnification factor associated with the root we have failed to find will be small. 
For every initial estimate, we obtain the roots i>(p). However, these roots do not necessarily 
differ; the different initial estimates may lead to an identical root. Therefore, in order to 
prevent multiple counting, we have to take paths that differ from each other. 



5.2. Determination of numerical parameters 

We determine appropriate numerical values of n and j/max to find all solutions v^p) of 
equation (!52|) whose magnification factors |/i(p)| significantly contribute to |-Ftotai|- If we fail 
to find such solutions, we underestimate the total fiux that emanates from the source to the 
observer. Hence, this is a very important task. Since the interference effect is not important 
for our present purpose, we neglect it for simplicity. Then the total flux is proportional to 
the magniflcation /i, deflned by 

p 
If we increase the values of n and ^max, the magniflcation fi will approach the maximum 
value /imax- We should determine the values of n and ymax so that the magniflcation /x 
satisfles fi ~ /Xmax- Since the lens distribution is determined using pseudorandom numbers, 
the magniflcation /i necessarily has a stochastic nature, and therefore we have to consider 
the ensemble average of fi. 

First, we set the source redshift zs to be 2 and numerically generate samples of lens 
distributions with numbers much larger than 20,000. Then, with the parameters |/max and 
n flxed, we calculate lens mappings to obtain the magniflcation fi for each lens distribution 
and construct an ensemble composed of 20,000 samples of fi. This procedure is repeated to 
construct ensembles for various values of parameters |/max and n. 
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We denote one of the ensembles by £, and further denote the ensemble average of a 
quantity Q{ii) related the magnification /i in £^ by {Q) g. The reader might think that (yu)^ 
will approach /imax in the limit of large |/max and n, and thus we might be able to choose 
the parameters T/max and n so that (/i)^ sufficiently converges to some value that will be 
equal to //max- However this is not the case, for the following reason. As shown in §5.3[ our 
numerical simulations imply that the probability P{^)d^ that the received magnification is 
in the range (/z, /i + djji) is proportional to /x~^ for sufficiently large /i. Therefore the second 
moment of the magnification (/i^), defined by 

/■oo 
(^2) ^ / ^^p[y,)dy, (55) 

will diverge. This implies that the law of large numbers fails, or in other words, (/i)^ cannot 
approximate the true avera ge value (u) even if the number of samples in the ensemble £ 



is much larger than unity (jPaczyhskil Il986l : iHolz fc Waldl Il998l : lRauchlll99ll ). Thus (/x) 



Is 



cannot be a criterion for determining whether the values of the parameters |/max and n are 
valid. 

To obtain the criterion for determining validity of the choices of ?/max and n, here we set 
an upper threshold /Xu on /i, and replace the magnification of samples in £ by /i^, if those are 
larger than /i^. We denote the ensemble obtained by this procedure as £. The probability 
distribution of the samples in the ensemble £ will be well approximated by 

/■oo 

P{^) = 9{fi^i - fi)P{fi) + S{fi - /iu) / P{ji)dji (56) 

if Umax and n are large enough, where 6{z) is a Heaviside's step function, and 5{z) is a Dirac 
delta function. Then we expect 

(/i)^ ~ / fiP{fi)dfi. (57) 

Jo 

We use (/i)^ as a criterion for the validity of the ensemble £. 

The relation between (/i)^ and the numerical parameters |/max and n is depicted in 
Figures [10] and [TTl where we have set /iu = 100. The value of (/i)g almost converges to the 
maximum values at ymax = 5,n = 500. In the case where the numerical parameters j/max 
and n are not sufficiently large, there are lens distributions {zi,yi) in which we cannot find 
any image. We exclude such cases as errors. In Table [1] we list the number of such errors 
appearing as we collected 20, 000 samples. Taking these results into account, hereafter we 
set ?/max = 5 and n = 500. In our results, the value of (/i)^ is smaller than the magnification 
in the global universe ^p. Here it should be noted that (/i)^ is not an observable quantity in 
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n=100}\ 




Fig. 10. — Ensemble average of magnification factor (/i)^, depicted as a function of the 
numerical parameters n and |/max- We have set zs = 2 and chosen the EdS universe as a 
global universe. The dashed horizontal line represent the value of iip- 
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Fig. 11. — Same as in Fig. ^U\, but here the global geometry is given by the FL universe with 
^AO = 0.7, i^mo = 0.3. 
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Table 1 : Number of errors in which we cannot find any image of the calculation for Figures 



[To] and [TT] is given here. 





fiAo = o,r 


mO — 


1 


^AO 


= 0.7,1] 


mO — 0.3 


n = 


100 X yrnax 


2/max ^ 5 


n = 100 


X 2/max 


1/max "^ 


2/max 


error 


n 


error 


i/max 


error 


n error 


0.5 


12 


10 


11791 


0.5 


6 


10 8495 


1.0 





20 


1294 


1.0 





20 864 


1.5 





30 


290 


1.5 





30 168 


2.0 





40 


83 


2.0 





40 61 


2.5 





50 


34 


2.5 





50 22 


3.0 





100 


2 


3.0 





100 1 


3.5 





200 


1 


3.5 





200 


4.0 





300 





4.0 





300 


4.5 





400 





4.5 





400 


5.0 





500 





5.0 





500 



the gravitational wave observations, since we can observe the only amplitude of superposed 
waves, although it is observable in optical observations, e.g., the observation of Type la 
supernov ae. As mention e d in ^ several works have d i scussed the equality between (u ) 
and /XF (|Weinbergl[l976l : IeUIs et Zlll998l : IcTaudell boool : Ir^ boOll : kibble fc Liei] I2OO5I ) . 
However, in practical sense, (/x) is not an observable quantity and thus (/i)^ might be more 
important. We will discuss this issue in connection with the optical observation elsewhere. 



5.3. Numerical results 



The distance-redshift relation of 2, 000 samples is plotted in Figures [I2] and [131 The 
probability distribution functions Pa(|-^| ) in several source redshifts are depicted in Figures 
[13] and [151 The results do not depend on the value of X/Ml if it is much smaller than unity. 
In the case of Pa(l-^P), the gradient in the high-amplification tail is more gradual than |-F|~^, 
whereas it is approximately proportional to \F\^ in the low-amplification tail (see Figs. [T4l 
and [TBI) . For comparison, we show the distance-redshift relation for incoherent waves (see 
Figs. [Tniand[T7|) and the probability distribution function of the magnification Pm^jj) (see 
Figs. [T8l and [T9|) . The probability of receiving an extremely high magnification is very small 
for each ray. Therefore, if a ray is highly magnified, it comes from only one lens and the lens 
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is located near the line of sight. In this case, we find 

yU oc y^^ for y ^ 
from equation (fT3l) . Hence, using equation (128|1 . we find 

Pm(/i) OC p (yiij,), z)-r '^ /^~^- 
We can see the consistent results in Figures [18] and [191 



(5^ 



(59) 




Redshift 

Fig. 12. — Distance- redshift relation in the case of short wavelength A ^ Ml. The global 
geometry is given by the EdS universe. We have numerically generated 2, 000 samples but 
we cannot plot all of the data points in this figure, since some samples have distances too 
large to plot within the frame. The DR distance is depicted by the solid line and the angular 
diameter distance in the global universe is represented by the dashed line. 



The distance dispersion is much larger than for long wavelengths, A ^ M^. The magnifi- 
cation of the amplitude comes from both the focusing and the interference effects, whereas the 
demagnification comes only from the interference effect. In contrast to the long- wavelength 
case, very large magnifications due to the focusing effect are possible for A -C M^, because 
there is no significant diffraction. As already mentioned, the first term on the right-hand side 
of equation (l50ll includes the shear an d focusing effects o n the ray bundle, which can also be 
found in the case of incoherent waves (JHolz &: Waldlll998l ). and these effects necessarily cause 
magnification of the amplitude (see Figs. [T6] and [T71) . By contrast, since the second term 
on the right-hand side of equation f[50l) can take either sign, the amplitude can be strongly 



-26- 



fim0=0-3.%0=0-7 



1 

0.9 
0.8 
0.7 
0.6 



Dyer- Boeder ■ 
Friedmann ■ 






Redshift 



Fig. 13. — Same as in Fig. [12], but here the global geometry is given by the FL universe with 
^AO = 0.7, fimo = 0.3. 



n„,n=i,a.n=o 



z,=0.5 
z,=1.0 
Zs=1.5 
zs=2.0 




Log(|Fr 



Fig. 14. — Probability distribution functions -Pad-f"! ) in several source redshifts. The global 
geometry is given by the EdS universe. We have numerically generated 100, 000 samples for 
each zs- 
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Fig. 15. — Same as in Fig. [T2|, but here the global geometry is given by the FL universe with 
^AO = 0.7, fimo = 0.3. 



^mO=l.^AO=0 




Fig. 16. — Distance-redshift relation in the case of incoherent waves with short wavelength 
A ^ Ml. We have numerically generated 2000 samples, taking into account only the focusing 
effect generated by the first term on the right-hand side of eq. (l50l) . The DR distance 
is depicted by the solid line and the angular diameter distance in the global universe is 
represented by the dashed line. 
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Fig. 17. — Same as in Fig. [16], but here the global geometry is given by the FL universe with 
^AO = 0.7, fimo = 0.3. 
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- Probability distribution functions PmifJ') in several source redshifts. The global 
is given by the EdS universe. We have numerically generated 100, 000 samples for 
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Fig. 19. — Same as in Fig. [12], but here the global geometry is given by the FL universe with 
^AO = 0.7, fimo = 0.3. 



demagnified, i.e., |Ftotai| < 1- Thus very large demagnifications are also possible for short 
wavelengths, A <^ Ml- As the result, the distance determined by the gravitational waves 
can grow longer than the DR distance. This is a peculiar nature of coherent waves. 

Here we evaluate the distance dispersion that could be obtained by the data provided by 
gravitational waves if our universe is well described by the present clumpy mass distribution 
model. For this purpose, we have to take a detection limit for weak events into account. Let 
us consider the observation of IHz gravitational waves emitted by neutron star-neutron star 
binaries that have the characteristic amplitude he, given by 



24/ 



9.88 X 10"^^ 1 + zs 



\5/6 



d] 



lOGpc 



-1 



i (cos^ e 

cos 9 



1) 



(60) 



wher e 9 and di are the inclination angle and luminosity distance of the binary (JHawking fc Israeli 
19891 ). respectively, and we have set the chirp mass of the binary to 1.2M0. The luminosity 
distance determined by lensed gravitational waves is given by (1 + zs)'^Djjji{0,zs)/ l-^totail- 
Assuming the planed sensitivity of DECIGO (ISeto et al.ll200lf ). we set the detection limit as 



K > 10 



-23 



(61) 



Randomly choosing the inclination angle 9, we numerically generate 20, 000 samples of he 
for each source redshift and then select those that satisfy equation fl6T|) . We denote the 
number of detectable samples as A'det- Results are shown in Tables [2] and [3] and Figures 
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and [2TI The distance dispersion SDohs is much larger than the estimate ob tained by 
assumin g a weak g r avita tional field for density perturbation (JTakahashil |2006| ) . This is 
because iTakahashil (120061 ) assumed a linear approximation for the lensed waveform even 
in the short wavelength limit. Furthermore, he considered the density perturbation as the 
inhomogeneity, but in this paper we consider all the matter in the universe to be the point 
masses, which is a more extreme inhomogeneous model. 



Table 2: Average distance (-Dobs) and distance dispersion 6Dohs for each source redshift, for 
EdS universe. The global geometry is given by the EdS universe. We have calculated the 
amplification factors |-Ftotai| of 20, 000 samples for each redshift. 



zs 


A^dct 


SDohs [Hq 


SDohs/ (-Dobs) 


(-Dobs) [Hq ] 


Ddr[Hq 


Df[H^'] 


0.25 


9339 


1.16x10-2 


6.81% 


1.71x10-1 


1.71x10-1 


1.69x10-1 


0.5 


8863 


3.23x10-2 


12.8% 


2.53x10-1 


2.55x10-1 


2.45x10-1 


0.75 


8450 


5.01x10-2 


16.9% 


2.96x10-1 


3.01x10-1 


2.79x10-1 


1.0 


8066 


6.48x10-2 


20.3% 


3.18x10-1 


3.29x10-1 


2.93x10-1 


1.25 


7673 


7.56x10-2 


22.9% 


3.29x10-1 


3.47x10-1 


2.96x10-1 


1.5 


7221 


8.35x10-2 


25.2% 


3.31x10-1 


3.60x10-1 


2.94x10-1 


1.75 


6780 


8.88x10-2 


27.1% 


3.28x10-1 


3.68x10-1 


2.89x10-1 


2.0 


6113 


9.32x10-2 


29.0% 


3.22x10-1 


3.74x10-1 


2.82x10-1 



Table 3: Same as Table [2l but here the global geometry is given by the PL universe with 
^AO = 0.7, i^mo = 0.3. 



zs 


iVdet 


SDohs [Hq 


SDohsl (Dohs) 


(-Dobs) [Hq ] 


Ddr [Hq 


DfWo'] 


0.25 


9307 


7.84x10-3 


4.15% 


1.89x10-1 


1.89x10-1 


1.88x10-1 


0.5 


8711 


2.41x10-2 


8.07% 


2.98x10-1 


2.99x10-1 


2.94x10-1 


0.75 


8079 


4.23x10-2 


11.7% 


3.63x10-1 


3.67x10-1 


3.53x10-1 


1.0 


7492 


6.03x10-2 


15.0% 


4.02x10-1 


4.11x10-1 


3.86x10-1 


1.25 


6862 


7.29x10-2 


17.2% 


4.24x10-1 


4.41x10-1 


4.02x10-1 


1.5 


6135 


8.41x10-2 


19.5% 


4.32x10-1 


4.61x10-1 


4.08x10-1 


1.75 


5408 


9.16x10-2 


21.2% 


4.33x10-1 


4.76x10-1 


4.07x10-1 


2.0 


4696 


9.43x10-2 


22.3% 


4.22x10-1 


4.87x10-1 


4.03x10-1 



The distance dispersion given in our results is lar ger than the typical value of observa- 
tional uncertainties evaluated without lensing effects (ICrowder fc Cornishll2005l ). However, 
our evaluation presented in this section is very basic one. In order to apply this evaluation to 
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Fig. 20. — Average distance and distance dispersion at each redshift in the clumpy universe. 
The global universe is given by the EdS universe. 
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Fig. 21. — Same as in Fi^ 
^AO = 0.7, i^mo = 0.3. 



Redshift 

but here the global geometry is given by the FL universe with 
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the future observational data, we need more detailed information about other experimental 
uncertainties and degeneracies between the situation presented here and others. We leave 
this for future works. 



6. Conclusion and summary 

We have studied the distance-redshift relation obtainable through gravitational waves 
in a clumpy universe. We have assumed that all of the matter takes the form of uniformly 
distributed point masses with identical mass, M^, and that the global geometry is described 
by the FL universe model. The distance from the source to the observer in the clumpy 
universe differs from that in the global universe, even if the redshift of the source is identical, 
since gravitational waves are magnified or demagnified by the gravitational lensing effect. We 
have numerically calculated the magnification due to multiple lensing effects in two extreme 
cases. In the first case, the wavelength A of the gravitational waves is much longer than the 
Schwarzschild radius, 2Mi, of the point-mass lens. In the second case, the wavelength A is 
much shorter than 2Ml. 

In the case of long wavelengths, A ^ Ml, the magnification due to a single lens is always 
very small, but it is not negligible. The total magnification in the multiple lens system is 
determined by the accumulation of all of the magnification effects due to each lens system, 
and is not small. We found that the relation between the average distance and redshift 
depends on the ratio A/M^, and it approaches the distance-redshift relation of the global 
universe as \/Mi increases. Furthermore, we found that the distance dispersion also depends 
on the ratio A/M^,, and it approaches zero as X/M^ increases. The distance-redshift relation 
in the global universe is found even in the clumpy universe if the wavelength of gravitational 
waves is much longer than the Schwarzschild radius of the lens objects. 

In the case of short wavelengths, A ^ M^, the geometrical optics approximation is 
valid. In order to calculate the total magnification, we employed the multiple lens plane 
method. We found that the distance dispersion is larger than that in the case of long 
wavelength. The total magnification comes from the focusing and interference effects on the 
ray bundles which emanate from the source to the observer. In this case, diffraction effects 
are negligible and thus the large magnifications due to the focusing effects on the ray bundles 
are possible. Moreover, interference effects can either enhance the magnifications or cause 
demagnificat ions . 

These lensing effects cause a distance dispersion in the distance-redshift relations ob- 
tained from gravitational wave data, and we have shown its properties. In particular, we 
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have clarified the dependence of tfie distance dispersion on tlie wavelengtli and tlie mass 
of the lens object. Making use of the properties of this distance dispersion, we can probe 
the inhomogeneities of our universe. The distance dispersion is sensitive to the ratio Mj^/A, 
and therefore our result suggest that we might use it to gain information about the typical 
masses of clumps that cause gravitational lensing effects. In order to obtain detailed infor- 
mation a bout the typical masses of lens objects, we need a method to calculate the case of 



Mi ~ A flYamamotol 120031 ). 



Suppose that we observe gravitational waves from binaries of compact objects in a 
narrow wave band centered at A. If the large distance dispersion does not appear in the 
distance-redshift relation obtained by these gravitational waves, we can impose a restriction 
on the typical mass scale of the clumping matters of M^ <C A. On the other hand, if a large 
distance dispersion appears, we can say that the mass of clumps Ml is much larger than the 
wavelength of the gravitational waves. 

In order to obtain a statistically reliable dispersion, we must detect a large number 
of gravitational waves. DECIGO^s planned sensitivity is 2 x lO^^^Hz"^' ^ around IHz; this 
mean s that 10^ chirp signals of coalescing binary neutron stars should be detected in 1 



year (ISeto et al.ll200ll ). Therefore, it should be possible to obtain reliable distance dispersion 
measurements. For example, if we observe gravitational waves from binary neutron stars in 
the waveband centered at IHz, we can determine whether the typical mass scale Ml of the 
dark matter that takes the form of macroscopic compact objects is much smaller or much 
larger than lO^M©. 

We are grateful to H. Ishihara and colleagues in the astrophysics and gravity group of 
Osaka City University for helpful discussion and criticism. 
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